#
# Produce `heatmap.pdf` from 2010 module data and estimates...
#  Run using R version 4.1.2
#
library(tidyverse) # Run using version: 1.3.1
library(haven)     # Run using version: 2.4.1

# Read in required estimates and data...
item_dat <- read_dta("../module2010_item_parameters.dta") %>%
                mutate(names = str_replace(names, "\\.", "_")) 
est_dat <- read_dta("../module2010_plus_estimates.dta") %>% 
              mutate( type = case_when( w1 > 0.5 ~ "Downsian",
                                        w2 > 0.5 ~ "Conversian",
                                        w3 > 0.5 ~ "Inattentive",
                                        TRUE ~ "Mixed"),
                       moderate = x > Hmisc::wtd.quantile(x, probs=1/3, weights=w1) & x < Hmisc::wtd.quantile(x, probs=2/3, weights=w1)) %>% 
              mutate(cces2008_nafta = ifelse(cces2008_nafta==1, 0, 1))

# Build long version (respondent-item) of data
long_dat <- est_dat %>% 
  select(respondent, moderate, type, item_dat$names) %>%
  pivot_longer(-(respondent:type), names_to="issue", values_to="response")

# Calculate fraction of one observed for each issue item
pct_ones <- long_dat %>%
              group_by(issue) %>%
              summarize(pct_ones = mean(response, na.rm=TRUE),
                        .groups="drop") %>%
              mutate(rank_ones = rank(-pct_ones)) %>%
              arrange(rank_ones)

# Make data displayed in the heat map..
gdat <- long_dat %>%
  # Make respondent-issue pair data 
  left_join(long_dat, by=c("respondent", "moderate", "type")) %>%
  filter(type =="Conversian" | (type=="Downsian" & moderate)) %>%
  group_by(type, issue.x, issue.y) %>%
  summarize(r00 = sum(response.x==0 & response.y==0, na.rm=TRUE),
            r01 = sum(response.x==0 & response.y==1, na.rm=TRUE),
            r10 = sum(response.x==1 & response.y==0, na.rm=TRUE),
            r11 = sum(response.x==1 & response.y==1, na.rm=TRUE),
             .groups="drop") %>%
  mutate( Odds = (log(r01+1) - log(r10+1)) ) %>%
  left_join(pct_ones %>% transmute(issue = issue, rank_ones.x = rank_ones), by = c("issue.x"="issue")) %>%
  left_join(pct_ones %>% transmute(issue = issue, rank_ones.y = max(rank_ones)+1-rank_ones), by = c("issue.y"="issue"))  

# Make the map
pdf("heatmap.pdf",height=5,width=10)
ggplot(gdat %>% mutate(type = ifelse(type=="Downsian", "Downsian Moderate", "Conversian")),
        aes(x=rank_ones.x, 
            y=rank_ones.y, 
            fill = Odds)) + 
  geom_tile() + 
  ylab("Items ordered by support for the liberal alternative") +  
  xlab("Items ordered by support for the conservative alternative") +
  facet_wrap(~type) + 
  theme_minimal() + 
  theme(legend.position="bottom") +
  coord_equal() + 
  geom_abline(intercept=133, slope=-1) +
  scale_fill_distiller(direction=-1, na.value = "white", oob=scales::oob_squish, 
                       breaks=c(  -rev(log(c(1, 3))),
                                  log(c(3))),  lim=c(-1.2, 1.2), 
                       labels = c("< 1:3" , "1:1",
                                   "> 3:1"), 
                                  palette="RdBu")
dev.off()